Is the broad histogram random walk dynamics 

correct? 

Jian-Sheng Wang 
Department of Computational Science, 
National University of Singapore, 
Singapore 119260 

2 October 1998 

Abstract 

We show explicitly that the broad histogram single-spin-flip random 
walk dynamics does not give correct microcanonical average even in one di- 
mension. The dynamics violates detailed balance condition by an amount 
which decreases with system size. As a result, in distribution different 
configurations with the same energy can have different probabilities. We 
propose a modified dynamics which ensures detailed balance and the his- 
togram obtained from this dynamics is exactly flat. 

Two years ago, Oliveira et al [j] |[ [| |j proposed an interesting method 
for Monte Carlo simulation of statistical mechanical systems. The number of 
degenerate states n(E) is computed from equations relating n(E) to the number 
of type of potential moves of a given dynamics. The Monte Carlo computation 
is purely combinatorial (or geometrical) unrelated to thermodynamics. The 
temperature dependence enters only after simulation in the weighting formulas. 
It is indeed an efficient method in comparison with histogram ^jj, especially 
multiple histogram method |(| . 

In the following, we use the Ising model to illustrate the method, but the 
idea generalizes easily to other models. One first chooses a protocol of Monte 
Carlo moves, for example, the single-spin flips. For a given state a with energy 
E = E(a) (a without an index denotes the set of all spins), consider all possible 
single-spin flips. The single-spin flips change the current state into TV possible 
new states, with new energy E' = E(a') = E + AE taking a small set of values. 
We classify the N new states according to AE, and count the number of such 
moves N(<j, AE). We have J2ae -^( cr : AE) = N. Then one can show rigorously 
that 1JJ 

n(E)(N(a, AE)) E = n(E + AE)(N(*', ~AE)) E+AE (1) 



as long as a set of moves and their reverse moves are both allowed. The average 
is a microcanonical average, i.e., average over all configurations having a fixed 
energy with equal weight: 
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Equation (|l|) is the fundamental result of broad histogram method. Once the 
quantity (N(cr, AE)) e is computed by any means, the density of states can be 
solved from Eq. (|l|). Since there are more equations than variables n(E) in 
general, one can either discard part of the equations, solving n(E) iteratively, 
or use least square method. In any case, once n{E) is determined from the set 
of equations, the canonical average is computed by 



There is no controversy on this part of the theory. However, the specific 
dynamics proposed in ref. JjJ is problematic. The dynamics is a random walk 
in energy as follows: pick a site at random, and compute energy change AE = 
E' — E if the spin is flipped. If the new state has a lower energy, E' < E, flip 
the spin. If the new state has a higher energy, AE > 0, flip the spin with a 
probability min(l, N(a, -AE)/N(a, AE)). If the energy is the same, AE = 0, 
flip the spin with some fixed probability, say 1. 

We believe that such a dynamics is wrong in the sense that it does not give 
correct microcanonical average with the samples generated in this dynamics. 
The stationary probability distribution of the dynamics depends not only on 
energy E but also on the partition N(a, AE) as pointed out already by Berg and 
Hansmann [Q . We show that such a dynamics does not satisfy detailed balance 
conditions with respect to an unknown but energy-dependent-only probability 
distribution. 

Let us first quote a mathematical theorem S| in the theory of discrete-time 
Markov chain using notations familiar to physicists. As before the symbol a 
will be a complete specification of the state (e.g., in Ising model, the set of all 
spins). Let Po(a) be the probability distribution of the initial configurations. If 
we start from a unique configuration (for example, from a ground state), Pq{&) 
will be 1 for the particular state, and for all other states. Let W(<t'|<t) be 
the transition probability, i.e., the probability that the system will be in state 
a' given that the current state is a\ W(a'\a) > 0, ^ , W(a'\a) = 1. Then the 
probability distribution after k steps of moves is 



Y,n(E)(A(a)) E exp(~E/kT) 



E 



(3) 



T 




E 




(4) 



2 



where W n = W ■ W ■ . . . W (n times) , viewing W as a matrix, and W n is matrix 
power. The theorem says that a unique large-step limit distribution 

P{&) = lim P n (a) = lim W n (a\a') (5) 

n — >oo n — >oo 

exists independent of the initial probability distribution, if 

1. W is irreducible, i.e., W k (a'\a) > for all states a and a' and some k > 0. 
In words, for any initial state a and any final state a', it is possible with 
nonzero probability to move from a to a' in k steps, k arbitrary (0, 1, 
2, ...) but finite. This is usually referred as ergodicity in some physics 
literature. 

2. W is aperiodic, i.e., W k (a\a) > for all k > k mini for all a. In words, 
starting from state a, the probability that it comes back to state a in k 
steps is nonzero for all k greater than a threshold value k m i n ■ This point is 
seldom mentioned in physics textbooks. When this condition fails, one can 
have probability distributions that oscillate between two or more forms. 
There will be no unique P(cr). 

Conditions 1 and 2 can be combined into a single statement: W k (a'\<r) > for 
all a and a' and for all k for sufficiently large k. 

The theorem guarantees that an equilibrium limit exists and is unique. If 
our aim is to simulate a known distribution, e.g., Boltzmann distribution, we 
also require that P{<j) is that distribution. It is sufficient then, in addition to 1 
and 2, W also satisfies 3 — detailed balance, 

W(<j\a')P(a') = W(a'\a)P(a), for all a and a'. (6) 

Example 1: The single-spin- flip Glauber rate with random selection of sites 
satisfies 1, 2, and 3 for temperature T > (but not if T = 0). In this case 

{0, if c and a' are not related by a single spin flip, 
^1 — tanh ^ , a' is obtained by a flip from a. ^ 

The diagonal elements W(ct|ct) are fixed by normalization (sum of the column 
of matrix W is 1, because probability of going anywhere is 1). 

Clearly condition 1 is satisfied — for arbitrary state a and a', we can move 
from a to a 1 by flipping one by one those sites which differ, we can do this in k 
steps where k is the number of spins which differ, condition 2 is also satisfied. 
Clearly VF(ct|(t) is nonzero for all a, so is the diagonal elements of all power of W. 
Condition 3 is satisfied with respect to the Boltzmann equilibrium distribution. 
Example 2: Consider a system of one spin. In each step, we flip the spin with 
probability one. Then 

W=(\ J), W* k =(l J), W 2k ^=W. (8) 
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This stochastic matrix violates condition 2. Thus, it docs not have a unique 
limit distribution. 

The broad histogram Monte Carlo dynamics in the basic form (single-spin- 
flip protocol) is as follows: 



W(a'\a) 



0, if a ^ Fa' and a ^ a'; 
±, if AS < and a = Fa'; ^ 

^min (l, N fa£g ), if AE > and a = Fa 



F changes configuration into a new configuration by a single-spin flip, and AE = 
E{a') — E(a) is the energy increment. VF(ct|ct) is again fixed by normalization 
(E<x' W(a'\a) = 1). The factor (1/N) on all of the terms reflects the fact that 
a site is picking at random. N is the total number of spins. 

The dynamics has problem with respect to condition 1, as the ground states 
are absorbing states. We cannot go to high energy states once the system is in 
the ground state. This small problem can be fixed by an ad hoc rule, e.g., if 
N(a, AE) = 0, AE < 0, reset it to infinity. 

Now, we look at condition 3, the detailed balance. Since we do not know 
P(a) for the above dynamics, we assume it depends on the state only through E, 
i.e., P(a) = f(E(a)) for some unknown function f(E). Then consider two ar- 
bitrary states a and a' related by a single-spin flip with energy E and E' , AE = 
E'-E> 0. From a to a', energy increases, W(a'\a) = (1/N) min(l, N(a, —AE)/ 
N(a,AE)). From a' to a, energy decreases, W(a\a') = 1/N. Now is 

f(E)(l/N)mm(l,N(a,-AE)/N(a,AE)) = f(E')(l/N) (10) 

for some function f(E)l Clearly, this is not possible in general unless the ratio 
N(a, —AE)/N(a, AE) is a function of E and E' only, and is independent of the 
current state a. Within the single-spin-flip protocol, the ratio does depend on 
a explicitly. We can easily find examples. 

P(a) can be computed in principle (and in practice for small systems) by 
solving WP = P, when the dynamics is precisely fixed. Since the configura- 
tions in a simulation appear with probability P(a), the average is ultimately 
related to average with the weight of P(a). Microcanonical average is restored 
as long as we have P(a) = f(E(a)) , i.e., the probability should be the same for 
configurations with the same energy. 

We present results of a simulation of a 5-spin chain with periodic boundary 
condition with the specific dynamics as discussed above. The data are presented 
in Table 1 together with exact results expected from true microcanonical aver- 
age. Monte Carlo steps (MCS) is 10 8 in the average, so we expect 4 digits of 
accuracy. 

The cases for E/J — —5 and 3 agree. But deviations of order 5% occur 
for the case E/J = —1. This is precisely due to the anomaly in P(a) in the 
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dynamics, since we compute (N(a, AE))e by an arithmetic average: 



(N(a,AE)) E = £ N(a,AE)P(a)/ ]T P(a') (11) 

ct, such that E(a)=E er\ E(<j')=E 

energy E, AE) /No. of samples. 

If P(cr) is a constant for a given E, then P(a) has no effect. But if P(a) for 
-E(er) = E, depends on a explicitly, then the above average is not microcanonical 
average. 

The specific unnormalized P(cr)s of the five-spin simulation dynamics are 
(computed with Mathematica exactly, which also agree with the simulation), 
P(+ + + + +) = 1, P(- + + + +) = 1, P(— +++) = 7/8, P(+- + -+) = 5/8. 
Other degenerate state probabilities are obtained by translation and spin up- 
down symmetry. We also have N(- + + + +AJ) = 2, N( h + +, 4 J) = 1, 

and N(- + + + +, -4 J) = 1, N( h + +, —4 J) = 0. From these numbers, 

we get (averages at £7/ J = — 1) 

W40> _LJ£^Wa = 1^.5333, (12, 

W-^))- '- 1 ^/" -!.^ (13, 

which is exactly what we get from Monte Carlo average and it is wrong! This 
feature is generic. Six and seven spin chains have the same problem. 

In actual implementation, a cumulative average of N(a, AE) is used instead 
of the instantaneous value for the transition rates 10. We have objections to 
the procedure. (1) If one uses the accumulated average (N(. . .)) up to the 
current state, then the dynamics is no longer Markovian — the transition rates 
not only depend on the current state but also on all the previous states. Then 
the theory of Markov chains does not apply. (2) Even if one uses the exact 
microcanonical average for the flip rate, the dynamics is still not correct in 
general. In Table 2 we present Monte Carlo computations of (N(cr, AE)) on a 4x 
4 lattice for (a) the original pure random walk dynamics, (b) the dynamics with 
the instantaneous value N(a, AE) replaced by exact microcanonical average 
values for the transition rates, and (c) the true microcanonical average by exact 
enumerations. Clearly, both the original dynamics and modified one do not give 
correct answers. 

In a variation of the method [Q , the dynamics is applied to an ensemble with 
energy limited to small window. By the same argument above, we see that this 
is also not correct. 

Detailed balance is a delicate matter. Since the detailed balance gives us 
more equations than the number of unknowns (P(<r)) in general, such equations 
are inconsistent unless the transition probabilities have special property. To 
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illustrate this point, let us take 



W(a'\a) = 6 a , Fa ,r(E{*')\E(a)) + d a , a ,r(E(a)\E{a)) , (14) 

where 5 a ,F<r' is 1 if er and a 1 differ by a single spin, and otherwise. This 
transition matrix corresponds to a single-spin-flip protocol with a transition 
rate r(E'\E) which depends on energies as an arbitrary function. Does this 
transition matrix preserve microcanonical property [-P(c) = /(-E(<r))]? The 
answer is not necessarily, unless 

r(E\E')f(E') = r(E'\E)f(E). (15) 
This can be shown as follows: 

P2(a) = ^WHOiV) (16) 
= r(E(a)\E(a))f(E(a)) + £iV(<7,i5' ' -E(a))r(E(a)\E') f(E') 

E' 

= f{E{a)) + 

J2 N(a, E'-E(a)) [r(E(a)\E')f(E') - r(E'\E(a))f{E(a)) . 

E> 

The last step used the fact that sum of the column elements in W is one. 
The first term depends on a through E implicitly as required. But there is 
no guarantee that the rest of the sum is a function of E only. However, when 
Eq. (|l5|) is satisfied, the rest of the terms are zero. 



We can eliminate the unknown function f(E) in Eq. (15) to get conditions 
on the matrix elements themselves. For example, for three distinct energies E, 
E', and E" with nonzero transition probabilities among them, we must have 

10, 

r{E\E")r{E"\E')r(E'\E) = r{E\E')r{E'\E")r{E"\E) > 0. (17) 

If there are four different energies with nonzero transition probabilities among 
them, we should have equations involving products of four of the transition 
probabilities. This would be the case in three-dimensional Ising model with 
single-spin flips. This equation is not satisfied in the random walk dynamics 
with i?-dependent-only rates where different size of energy jumps are allowed. 

From the results presented in refs. Q [|, ||, || the deviations for thermody- 
namics quantities for large systems are almost unnoticeable. Let us discussion 
in the context of Eq. jl7| ) why this is so. For the two-dimensional Ising model, 
this equation is the only relevant constraint equation. We define 



v{E) 



{ r(E\E")r(E"\E')r(E'\E) 



r{E\E')r{E'\E")r{E"\E) 



(18) 



G 



as the detailed balance violation, where we take E' = E + 4J, and E" = E + 8J. 
For the random walk dynamics, 

r(E ' |E, = {^(i.|Sfefc*)- (19) 

Substituting this expression into Eq. (18), we obtain (for E < 0) 
(N{a,-U)) E (N(a',-U)) Et (N(a,8J)) E 



v(E) = 



1 - 



(N(a, U)) E (N(a>, AJ)) B > (N(a, -8J))j 



(20) 



If we take the large-size limit, then the functions (N(- ■ •)) are smooth functions 
in energy, and we approximate the discrete spectrum by continuous functions, 
rii{u) = limjv_ >00 (l/AT)(Ar(o-, iAJ)) uN , i — 0, ±1,±2, u = E/N. It was pointed 
out to us by Oliveira jO that the functions n;(u) can be related to thermody- 
namics, sec also rcf. |l2|| . Here we give a slightly different argument for it. Since 
rii(u) are large-size microcanonical average, using the equivalence between dif- 
ferent ensembles in the thermodynamic limit, we can compute the same qualities 
by canonical ensemble. rii(u) is simply the probability that a site is surrounded 
by i + 2 spins of the same signs. Evoking Boltzmann distribution, we can show 
rii(u)/n-i(u) = exp(i4J/3(it)) , where /3(u) = 1/kT, T is the temperature for 
system at an average energy u per spin. Using this result, Eq. (|2C| ) can be 
further simplified [|T| netly as 



v(u) 



1 _ e ij(f3(u)-l3(u+AJ/N)) 



dp 16J 



2 



du N 



(21) 



This result suggests that the violation is of order 1/N in general. At the 
critical region, it becomes even better, by a factor of the inverse specific heat. 
For the two-dimensional Ising model at T c (near u — u c — — V%), we have 
v(u c ) oc l/(N\ogN). This indeed justifies the results obtained by computer 
simulation on large system. If instantaneous value N(a, AE) is used, we expect 
errors of order 1 / VN. 

Since detailed balance violation is a small perturbation to the transition rate, 
we expect the Monte Carlo result for (N(- ■ •)) also violates detailed balance by 
an amount proportional to v(u). We check this, replacing r(E'\E) by Monte 
Carlo estimates of (N(cr,E' — E))e in Eq. (|l8|). Let us call this quantity v(u). 
Numerically we found v(u c ) = 0.3,0.12,0.03,0.008 for two-dimensional square 
lattice systems of linear size L = 4, 8, 16, 32, respectively. The results are 
consistent with Eq. (^l|). 

However, we do not like approximate algorithms. Finite-size scalings are 
important part of computer simulation methods. It is thus desirable to have 
exact algorithms with any lattice sizes. A remedy to the problem is to use 
a single function f(E), similar to ref. |l3| (which uses f(E) on l/n(E)). For 
example, let f(E) = l/(N(a,AE)) E for a fixed AE. Use a flip rate r = 
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min[l, f(E(a'))/f(E(a))] for a move from a to a'. This will at least fulfill the 
detailed balance. The second possibility is to restrict energy changes to only 
three possible values 0, ±4J. The flip rates can be any arbitrary function of 
E and AE = 0, ±4 J. When such a condition is imposed, we can always find a 
consistent solution f(E) to the detailed balance equations. This will introduce 
ergodicity problem, but inconsistency to equations like (|l7|) will not occur. This 
also explains why the dynamics works in one dimension if iS-dependent-only flip 
rates are used. 

The last and perhaps also the best recommendation is to take 

without restriction to energy changes. When this rate is used, Eq. ( p7[ ) and simi- 
lar constraint equations are satisfied automatically (as long as the quantities are 
exact microcanonical averages). This is so because (N(a, AE))e themselves are 



also transitions rates [|8| and also satisfy equations similar to Eq. (17). Moreover, 
we known exactly what is the stationary distribution f(E). It is just the inverse 
density of states, .i.e., f(E) oc l/n(E). The detailed balance conditions for this 
new dynamics are simply equivalent to Eq. (|f|). The dynamics also gives an 
exactly flat histogram. The histogram in general is proportional to n(E)f(E). 
Since f(E) oc l/n(E), the energy histogram is a constant. 

For actual implementation, we must determine the microcanonical average 
iteratively, using at first cumulative average of N(<j, AE) for transition rates. 
After a first approximation, we can either refine them further, or use the slightly 



incorrect one but reinforce Eq. (17) exactly on data. When this is done, the 
histogram will not be exactly flat, but the dynamics is correct. Detail study of 
this dynamics will be presented elsewhere. 

In conclusion, we have shown that the original random walk dynamics is only 
an approximation. Although it gives good results for large systems, on small 
systems, the value for (N(- ■ •)) can have errors as large as 20%. We proposed new 
transition rates that satisfy detailed balance condition. In particular, instead 
of using the ratio of number of moves to states with energy ±AE from the 
current state, we use the ratio of average number of moves from new states to 
the current energy states and from current state to the new states. This very 
small change to the original method plus the requirement that the average values 
should be used ensures detailed balance. In addition, the histogram obtained 
by this dynamics is exactly flat. 

The author thanks R. H. Swcndscn for drawing attention to this problem. He 
particularly thanks P. M. C. de Oliveira for pointing out errors in the derivations 
of Eq. (^o|) and (^l|), for many clarifications, arguments, stimulating discussions, 
and critical readings of several versions of the manuscripts. 
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Table 1: The average value (N(a, AE)) for a system of five-spin chain in broad 
histogram random walk dynamics from Monte Carlo simulation and exact results 
with true microcanonical average. 



E/J 


(N(a,AE = -AJ)) E (N(a, AE = +AJ)) E 


Random Walk Dynamics Result 


-5 


5 


-1 


0.53326 1.53326 


3 


3 


Exact Microcanonical Results 


-5 


5 


-1 


1/2 3/2 


3 


3 



Table 2: Monte Carlo random walk dynamics and exact results on a 4 x 4 square 
lattice with periodic boundary conditions: (a) pure broad histogram random 
walk dynamics, using instantaneous N(a,AE), MCS = 1.45 x 10 9 ; (b) broad 
histogram random walk dynamics with strictly independent flip rates, using 
the exact (N(a,AE)) E , MCS = 1.45 x 10 10 ; (c) true microcanonical average 
computed by exact numeration of all the states. Statistical/numerical error is 
on the last digit. 



E/{AJ) 


(N(a,AE=-8J)) E 


(N{a,AE=^J)) E 


(N{a,AE = 0)) E 


(a) Random walk dynamics, current N(a, AE) 


4 


0.9878 


0.7850 


1.6194 


5 


0.4348 


1.5815 


3.0233 


6 


0.4077 


1.7295 


5.4252 


7 


0.5626 


2.6570 


5.9934 




(b) Random walk using average (N(a, 


AE)) 


4 


0.86573 


0.80220 


1.6888 


5 


0.39313 


1.5878 


3.2620 


6 


0.47274 


1.7930 


5.3393 


7 


0.48124 


2.9319 


5.8429 


(c) Exact microcanonical average 


4 


0.83018868 


0.90566038 


1.81132075 


5 


0.29629630 


1.55555556 


3.55555556 


6 


0.38755981 


1.81818182 


5.51196172 


7 


0.45283019 


2.94339623 


5.88679245 
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